Crystallization of a classical two-dimensional electron system: 
Positional and orientational orders 
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Crystallization of a classical two-dimensional one-component plasma (electrons interacting with the 
Coulomb repulsion in a uniform neutralizing positive background) is investigated with a molecular- 
dynamics simulation. The positional and the orientational correlation functions are calculated, to 
the best of our knowledge, for the first time. We have found an indication that the solid phase 
has a quasi-long-range (power-law) positional order along with a long-range orientational order. 
This indicates that, although the long-range Coulomb interaction is outside the scope of Mermin's 
theorem, the absence of ordinary crystalline order at finite temperatures applies to the electron 
system as well. The 'hexatic' phase, which is predicted between the liquid and the solid phases by 
the Kosterlitz-Thouless-Halperin-Nelson- Young theory, is also discussed. 
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Wigner pointed out in the 1930's that an electron 
system should crystallize dujC to the Coulomb repul- 
sion for low enough densities .EJ Although quantum effects 
are essential, the concept of electron crystallization can 
be generalized to classical cases. In fact, Grimes and 
Adams succeeded in observing a liquid-to-solid transi- 
tion in a classical two-dimensional|i2D) electron system 
on a liquid-helium surface in 1979.13 In this system elec- 
trons obey classical statistics because the Fermi energy 
is much smaller than ksT. The thermodynamic prop- 
erties of the classical electron system are wholly deter- 
mined by the dimensionless coupling constant F, the ra- 
tio of the Coulomb energy and the kinetic energy. Here 
F = {e^/4:TTea)/kBT, where a ~ (vrn)"^/^ is the mean- 
distance between the electrons, n the density of electrons, 
e the charge of an electron, and e is the dielectric con- 
stant of the substrate. For F ^ 1 the system will behavgj 
as a gas while for F ^ 1 as a solid. Grimes and AdamsBI 
found that the phase transition occurs at Fc = 137 ± 15. 

On the other hand, Mermin proved rigorously, more 
than thirty years ago, that no true long-range crystalline 
orders are possible ia the thermodynamic limit at finite 
temperatures in 2D .13 In the proof, short-range interac- 
tions are assumed, while the 1/r Coulomb interaction is 
too long-ranged to apply Mermin's arguments (see Ref. g 
for more precise mathematical conditions!- Although 
there have been some theoretical attemptafl'El to extend 
the theorem to the long-range Coulomb case, no rigorous 
proof is attained up to now. Numerically, Gann et al.tl 
investigated the problem with a Monte Carlo method by 
calculating the root- mean-square displacement. However 
they found it difhcult to rule out the possibility that 
the root-mean-square displacement approaches a con- 
stant value in the thermodynamic limit, which is a con- 
ventional definition of a solid. 

Another intriguing problem concerning 2D systems is 
the melting mechanism, where the Kosterlitz-Thouless- 
Halperin-Nelson- Young (KTHNY) theoryD'El has pre- 
dicted the existence of the 'hexatic' phase between the 
liquid and the solid phases (see Refs. H, |l^, ^ for re- 



views). The hexatic phase is characterized by a short- 
range positional order and a quasi-long-range orienta- 
tional order. Numerical calculations have been carried 
out by spjexaLauthors with a molecular-dynamics (MD) 
method.OEjO Most of the authors have obtained the 
melting point in good agreement with the Grimes- Adams 
experiment. However the results are rather contcover- 
sial on the verification of the hexatic phase. Morili3 con- 
cluded, from a MD result for the shear modulus, that the 
result agrees well witii the KTHNY theory. On the other 
hand, Kalia et al.t3 observed hysteresis in the tempera- 
ture dependence of the total energy in their MD simula- 
tion to conclude that the melting is a first-order phase 
transition, which is incompatible with the KTHNY the- 
ory. The Monte Carlo study by Gann et alfl is also in- 
dicative of a first-order phase transition. 

In addressing both of the above problems, i.e., the 
range of the ordering in the solid phase and the nature of 
the melting, the most direct way is to investigate the po- 
sitional and the orientational correlation functions, since 
the phases should be characterized in terms of them. This 
is exactly the purpose of the present paper, which is done, 
to the best of our knowledge, for this system for the first 
time. 

In order to accurately incorporate the temperatur| 
have employed Nose-Hoover's canonical MD methoc 
for the classical 2D electron system, while the previ- 
ouSpCalCjUlations were done for micro-canonical ensem- 
blc.EJOEj Electrons are confined to a rectangle with a 
rigid uniform neutralizing positive background. We im- 
pose periodic boundary conditions. An electron then in- 
teracts with infinite arrays of the periodic images with 
the long-range l/r potential. We employ the Ewald sum- 
mation method to take care of this. The rectangle is 
chosen to be close to a square to minimize surface ef- 
fects. The aspect ratio of the rectangle is taken to be 
Ly/Lx — 2/>/3, which can accommodate a perfect trian- 
gular latticetZl with N = 4M^ (M: an integer) particles. 
The equations of motion are integrated numerically with 
Gear's predictor-corrector algorithm. We adopt a time 



step of 1.0 X 10~^^ sec, which guarantees six-digit ac- 
curacy in the energy conservation after several tens of 
thousands of steps. 

We have performed the simulation both from a typi- 
cal liquid (r — 60) and from a typical solid (F — 200). 
The initial conditions are set as follows: The electrons 
are placed randomly in the liquid phase or placed at the 
perfect triangular lattice points in the solid phase. The 
velocities of the electrons are assigned according to the 
Maxwell-Boltzmann distribution in either case. The low- 
est energy configurations are sought with a simulated an- 
nealing method. Namely, the positions and the velocities 
of the electrons are updated for a certain time interval. 
Once a thermal equilibrium sets in, the temperature is 
raised or lowered by a small amount. The latest po- 
sitions and velocities are used as the initial conditions 
for the simulation at the new temperature. This proce- 
dure temporarily puts the system out of equilibrium, but 
updating the positions and velocities for a certain time 
interval equilibrates the system. This annealing process 
is repeated from a liquid phase to a solid phase (or vice 
versa) across the transition. We take care that the sys- 
tem is well equilibrated, especially near and after the 
transition, by allowing large numbers of time steps. The 
results presented in this paper are for N = 900 electrons 
with MD runs from 30,000 to 110,000 time steps for each 
value of r. The correlation functions are calculated for 
the last 20,000 time steps, j— , 

Following Cha and Fertigjl3 we define the positional 
and the orientational correlation functions from which 
we identify the order in each phase. First, the positional 
correlation function is defined by 

C(r) ^ (p^(r)pG(O)) (1) 

V (2) 



^5(r - |r, - V.J 



where G is the reciprocal vector of the triangular lat- 
tice, and PG(r) — exp(iG • r). The angular brackets in 
Eq. (|I|) stand for both the summation over particles and 
the thermal average. In Eq. (H) a summation is taken 
over six reciprocal vector G's that give the first peaks of 
the structure factor. In practice, the ^-function must be 
broadened so that it can be handled numerically. 
The orientational correlation function is defined by 
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where V'eCr) — -^ X^a" e^'^°(''\ and 6*0, (r) is the angle of 
the vector connecting an electron at r and the a-th near- 
est neighbor with respect to a fixed axis, say, the a;-axis. 



The summation is taken over Uc nearest-neighbors which 
are determined by the Voronoi diagrarr£3, or equivalcntly 
its dual mapping, the Delaunay triangulation. 

We first look at the positional and the orientational 
correlation functions for F = 200 and F — 160, typical 
solid phases, in Fig. |l|. Although F (T) is high (low) 
enough, the positional correlation is seen to decay slowly 
in both cases, indicating an algebraic decay at large dis- 
tances. The round-off in the correlation function around 
half of the linear dimension of the system size is consid- 
ered to be an effect of the periodic boundary conditions. 
The algebraic decay of the positional correlation function 
implies that the 2D electron solid has only &jmasi-long- 
range positional order at finite temperatures .B2l Thus we 
have [Obtained a numerical indication that Mermin's the- 
oremEl applies to the electron system as well, which is 
consistent with the analytical (but not rigorous) results 
obtained in Refs. H, pi 

By contrast, the orientational order is seen to be long- 
ranged. Therefore while the 2D electron solid has no true 
long-range crystalline order, it is topologically ordered. 
The triangular structure is seen as the peaks in both the 
positional and the orientational correlation functions (see 
the inset of Fig. IT]) . 

We have plotted the Delaunay triangulation of a snap- 
shot of the electron configuration for a solid (F — 160) or 
for a liquid (F — 90) in Fig. 0. We can in particular look 
at the topological defects, i.e., five-fold and seven- fold co- 
ordinated electrons. From the result the defects are seen 
to appear in isolated pairs, or more precisely in quartets, 
in the solid phase, which explains how a quasi-long-range 
positional order is compatible with a long-range orienta- 
tional one. On the other hand, defects appear with a 
high density in the liquid phase. 

We now focus on the orientational correlation function 
near the crystallization in Fig. ||. The result is obtained 
by cooling the system from a liquid to a solid. The ori- 
entational order is short-ranged for F < 120 and long- 
ranged for F > 140. Around the liquid-solid boundary 
(F = 130), the orientational correlation function, plot- 
ted on a double logarithmic scale in Fig. |[ indicates 
an algebraic decay (while the positional order is short- 
ranged). The power evaluated from the data is approx- 
imately equal to unity, which is greater than the upper 
bound of 1/4 predicted by the KTHNY theory. Large 
statistical errors near the transition, however, prevent us 
from drawing any definite conclusion on the existence of 
the hexatic phase. In fact, the correlation functions be- 
have like a solid at F = 130 when the system is heated 
from a solid to a liquid with no indication of the hexatic 
phase. This may be due to a finite-size effect, where a 
solid can be pinned in a melting process. For a small 
system, N = 100, a solid phase in fact persists down to 
F = 120 when heated. 

Another finite-size effect is that, when the system crys- 
tallizes, the crystal axes can tilt from the unit cell axes 
of the finite system. The crystallization does occur in a 
tilted way in the present simulation. The misalignment 



causes a long relaxation time for the system to reach 
the lowest energy state. However the fact that the sys- 
tem crystallizes with tilted axes shows in itself that the 
N — 900 system is sufficiently large in that boundary 
effects are not too strong. By contrast, we found that 
the crystal axes are always aligned to the unit cell for 
N = 100. For the N =pLpO system, which is the size 
employed by Kalia et al.E3, we found no indication of 
the hexatic phase, either. A numerical difficulty in MD 
simulations also arises from finite time steps. Finite time 
effects might result in insufficient equilibration, especially 
near a continuous transition. Even if the system is well- 
equilibrated, it would be difficult to tell a slow expo- 
nential decay from an algebraic decay in the correlation 
function for a finite system. 

In summary, we have performed a molecular-dynamics 
simulation to investigate the ordering of a classical 2D 
electron system. From the positional and the orienta- 
tional correlation functions we have found an indication 
that there is a quasi-long-range positional order and a 
long-range orientational order in the solid phase, which 
implies that Mermin's theorem is not spoiled even for the 
long-range 1/r interaction. On the other hand, we have 
obtained only a sign, not a conclusive result, for the exis- 
tence of the hexatic phase predicted by KTHNY theory, 
which thus remains an open question. 

Although wc have basically electrons on a liquid- 
helium surface in mind, a planar classical one-component 
plasma is recently realized as laser-cooled ions trapped in 
a disk region.EJ The disk has a finite thickness, for which 
stable crystalline phases are observed. If the ions could 
be trapped completely in 2D, the present picture would 
be applicable. Conversely, it is an interesting theoretical 
problem to extend the present line of approach to planar 
systems with finite thicknesses. 

We wish to thank Kazuhiko Kuroki, Hiroshi Imamura, 
Katsunori Tagami, and Naruo Sasaki for valuable discus- 
sions. The numerical calculations were mainly done with 
Fujitsu VPP500 at the Supercomputer Center, Institute 
for Solid State Physics, University of Tokyo. 
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FIG. 1. The positional (the upper panel) and the orienta- 
tional (the lower panel) correlation functions for F = 200 and 
F — 160. The horizontal scale is in units of the lattice con- 
stant of the triangular lattice. The inset in the lower panel 
illustrates the first four distances which give the peaks in the 
correlation functions. 



FIG. 2. The Delaunay triangulation of a snapshot of the 
electron configuration for a solid with F — 160 (a) or for a 
liquid with F = 90 (b). Five(seven)-fold coordinated electrons 
are marked with empty (solid) circles. In (b) four(eight)-fold 
coordinated electrons are marked with empty (solid) squares. 
A part of the system is shown in either panel. 
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FIG. 3. The orientational correlation function near the 
crystallization. The result was obtained by cooling the system 
from a liquid to a solid. 
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